home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / Scientist's Helper src / s.helper.5 / regspline.c < prev    next >
C/C++ Source or Header  |  1986-02-06  |  8KB  |  317 lines
  1. #include "all.h"
  2. #include "regtabext.h"
  3.  
  4. LoadSg3() /*force load of segment*/
  5. {
  6.     ;
  7. }
  8.  
  9. Spline( X,Y, N, C, X0,Y0, F, AM,BM,CM,DM,EM,FM, LAST )  /* natural cubic splines */
  10. float X[], Y[], C[], X0, *Y0, AM[], BM[], CM[], DM[], EM[], FM[];
  11. int N, F, *LAST;
  12.  
  13. /*     X        THE ARRAY OF X VALUES ARRANGED IN ORDER OF INCREASING X VALUE WITH NO DUPLICATIONS
  14.     Y        THE ARRAY OF Y VALUES
  15.     N        THE NUMBER OF (X,Y) PAIRS
  16.     C        AN ARRAY INTO WHICH THE SPLINE COEFFICIENTS ARE STORED
  17.     X0        AN ARBITRARY X VALUE AT WHICH THE SPLINE IS TO BE EVALUATED
  18.     Y0        POINTER TO RETURN THE VALUE OF Y AT X0
  19.     F        A FLAG, COMPUTE COEFFICIENTS ON TRUE, EVALUATE SPLINE AT X0 ON FALSE
  20.     AM        TEMP ARRAY OF LENGTH N+2
  21.     BM        TEMP ARRAY OF LENGTH N+2
  22.     CM        TEMP ARRAY OF LENGTH N+2
  23.     DM        TEMP ARRAY OF LENGTH N+2
  24.     EM        TEMP ARRAY OF LENGTH N+2
  25.     FM        TEMP ARRAY OF LENGTH N+2
  26.     LAST    POINTER TO INTEGER FOR TEMPORARY STORAGE */
  27. {
  28.     int I, J, JP1;
  29.     float HJ;
  30.     char s[80];
  31.             
  32.     if (F) {  /* compute spline coefficients */
  33.         *LAST=1;
  34.         for (I=2; I<=(N-1); I++ ) {
  35.             X0 = (Y[I+1]-Y[I])/(X[I+1]-X[I]);
  36.             X0 = X0 - (Y[I]-Y[I-1])/(X[I]-X[I-1]);
  37.             DM[I-1] = 6.0*X0;
  38.             BM[I-1] = 2.0*(X[I+1]-X[I-1]);
  39.             AM[I-1] = X[I+1] - X[I];
  40.             }
  41.             
  42.         TRIDI(AM,BM,CM,DM,C,N-2,EM,FM);
  43.         
  44.         C[N]=0.0;  /* SHIFT ALL THE C'S AND SET C(1)=C(N)=0.0 */
  45.         for( I=2; I<=(N-1); I++ ) {
  46.             J = N-I;
  47.             C[J+1]=C[J];
  48.             }
  49.         C[1] = 0.0;
  50.         
  51.         } /*end if compute spline coefficients*/
  52.         
  53.         
  54.     else { /* evaluate spline at X0 */
  55.         
  56.         /* locate proper knot */
  57.         if( (X0<=X[(*LAST)+1]) && (X0>=X[*LAST]) ) { /*check current knot*/
  58.             J = *LAST;
  59.             JP1 = J+1;
  60.             } /*end if*/
  61.         else { /*check elsewhere*/
  62.             (*LAST)++;
  63.             if ( (*LAST)>(N+1) ) *LAST=N+1;
  64.             if ( (X0<=X[(*LAST)+1]) && (X0>=X[*LAST]) ) { /*check next knot*/
  65.                 J = *LAST;
  66.                 JP1 = J+1;
  67.                 }
  68.             else if (X0<X[1]) { /*check off beginning*/
  69.                 J = 1;
  70.                 JP1 = 2;
  71.                 }
  72.             else if (X0>X[N]) { /*check off end*/
  73.                 J = N-1;
  74.                 JP1 = N;
  75.                 }
  76.             else { /*scan all knots*/
  77.                 for( I=1; I<=(N-1); I++ ) {
  78.                     if( (X0<=X[I+1]) && (X0>=X[I]) ) {
  79.                         J=I;
  80.                         JP1=I+1;
  81.                         break;
  82.                         } /*end if in interval*/
  83.                     } /*end for*/
  84.                 } /*end else*/
  85.             } /*end if*/
  86.         *LAST = J;
  87.         
  88.         /* COMPUTE Y0 AT X0 */
  89.         HJ = X[JP1] - X[J];
  90.         *Y0 = (C[J]/(6.0*HJ)) * (X[JP1]-X0)*(X[JP1]-X0)*(X[JP1]-X0);
  91.         *Y0 += (C[JP1]/(6.0*HJ)) * (X0-X[J])*(X0-X[J])*(X0-X[J]);
  92.         *Y0 += (Y[JP1]/HJ-C[JP1]*HJ/6.0) * (X0-X[J]);
  93.         *Y0 += (Y[J]/HJ-C[J]*HJ/6.0) * (X[JP1]-X0);
  94.         } /*end if evaluate spline*/
  95. }
  96.  
  97. SplInterpolate()
  98. {
  99.     float *xd, *yd, *c, *am, *bm, *cm, *dm, *em, *fm; /*temporary arrays*/
  100.     float firstX, lastX, x0, y0;
  101.     int npair, col, i, j, k, newRows, oldRows, last, goodX;
  102.     long bytesNeeded;
  103.     char str[cmdWordLen];
  104.  
  105.     if( table.header.interpolated ) {
  106.         ErrMsg("table already interpolated");
  107.         }
  108.     if( table.header.rows < 4 ) {
  109.         ErrMsg("not enough rows");
  110.         }
  111.     if( table.header.cols < 2 ) {
  112.         ErrMsg("not enough cols");
  113.         }
  114.         
  115.     SToR( command.cmdWord[1], &(table.header.samp), FALSE );
  116.     if (table.header.samp<=0)  {
  117.         table.header.samp = 1.0;
  118.         ErrMsg("samp must be > 0");
  119.         }
  120.     
  121.     /* find first non-NaN x value */
  122.     j=FALSE;
  123.     for( i=1; i<=table.header.rows; i++ ) {
  124.         GetTable( i, 1, &firstX );
  125.         if( !NaN(&firstX) ) {
  126.             j=TRUE;
  127.             k=i;
  128.             break;
  129.             } /*end if*/
  130.         } /*end for*/
  131.     if( (!j) || (k==table.header.rows) ) ErrMsg("not enough data in col 1");
  132.     
  133.     /* find last non-NaN x value */
  134.     for( i=table.header.rows; i>=1; i-- ) {
  135.         GetTable( i, 1, &lastX );
  136.         if( !NaN(&lastX) ) {
  137.             break;
  138.             } /*end if*/
  139.         } /*end for*/
  140.         
  141.     /* count up non-NaN x values */
  142.     goodX=0;
  143.     for( i=1; i<=table.header.rows; i++ ) {
  144.         GetTable( i, 1, &x0 );
  145.         if( !NaN(&x0) ) {
  146.             goodX++;
  147.             } /*end if*/
  148.         } /*end for*/
  149.     if( goodX<4 ) ErrMsg("not enough data in col 1");
  150.  
  151.     /* check for points that are out of order */
  152.     x0 = firstX;
  153.     for( i=k+1; i<=table.header.rows; i++ ) {
  154.         GetTable( i, 1, &y0 );
  155.         if( !NaN(&y0) ) {
  156.             if( y0<=x0 ) ErrMsg("data in col 1 out of order");
  157.             x0 = y0;
  158.             } /*end if*/
  159.         } /*end for*/
  160.             
  161.     bytesNeeded = 4L * (4+goodX);
  162.     if( (xd=NewPtr(bytesNeeded)) == 0L ) {
  163.         ErrMsg("not enough free memory");
  164.         }
  165.     if( (yd=NewPtr(bytesNeeded)) == 0L ) {
  166.         DisposPtr(xd);
  167.         ErrMsg("not enough free memory");
  168.         }
  169.     if( (c=NewPtr(bytesNeeded)) == 0L ) {
  170.         DisposPtr(xd);
  171.         DisposPtr(yd);
  172.         ErrMsg("not enough free memory");
  173.         }
  174.     if( (am=NewPtr(bytesNeeded)) == 0L ) {
  175.         DisposPtr(xd);
  176.         DisposPtr(yd);
  177.         DisposPtr(c);
  178.         ErrMsg("not enough free memory");
  179.         }
  180.     if( (bm=NewPtr(bytesNeeded)) == 0L ) {
  181.         DisposPtr(xd);
  182.         DisposPtr(yd);
  183.         DisposPtr(c);
  184.         DisposPtr(am);
  185.         ErrMsg("not enough free memory");
  186.         }
  187.     cm = am-1;
  188.     if( (dm=NewPtr(bytesNeeded)) == 0L ) {
  189.         DisposPtr(xd);
  190.         DisposPtr(yd);
  191.         DisposPtr(c);
  192.         DisposPtr(am);
  193.         DisposPtr(bm);
  194.         ErrMsg("not enough free memory");
  195.         }
  196.     if( (em=NewPtr(bytesNeeded)) == 0L ) {
  197.         DisposPtr(xd);
  198.         DisposPtr(yd);
  199.         DisposPtr(c);
  200.         DisposPtr(am);
  201.         DisposPtr(bm);
  202.         DisposPtr(dm);
  203.         ErrMsg("not enough free memory");
  204.         }
  205.     if( (fm=NewPtr(bytesNeeded)) == 0L ) {
  206.         DisposPtr(xd);
  207.         DisposPtr(yd);
  208.         DisposPtr(c);
  209.         DisposPtr(am);
  210.         DisposPtr(bm);
  211.         DisposPtr(dm);
  212.         DisposPtr(em);
  213.         ErrMsg("not enough free memory");
  214.         }
  215.         
  216.     table.header.start = firstX;
  217.     x0 = 1.0 + ((lastX-firstX)/table.header.samp);
  218.     oldRows = table.header.rows;
  219.     if ( ((long)x0) > ((long)table.header.maxRows) ) {                   
  220.         WriteLine("warning. some data lost from end of columns");
  221.         newRows = table.header.maxRows;
  222.         } /*end if*/
  223.     else {
  224.         newRows= (int)x0;
  225.         } /*end if*/
  226.  
  227.     for (col=2;  col<=table.header.cols; col++ ) {
  228.     
  229.         /* stuff arrays xd and yd with non-Nan data */
  230.         i=1;
  231.         npair=0;
  232.         while( NextNotNan( i, 1, col, &i ) ) {
  233.             npair++;
  234.             GetTable( i, 1,   (xd+npair) );
  235.             GetTable( i, col, (yd+npair) );
  236.             i++;
  237.             } /*end while*/
  238.             
  239.         if (npair<4) {
  240.             IToS( col, str );
  241.             WritePhrase("warning: not enough data in col ");
  242.             WriteLine(str);
  243.             table.header.rows = newRows;
  244.             for( j=1; j<=table.header.rows; j++ ) {
  245.                 SetTable(j,col,infinity,FALSE);
  246.                 }
  247.             table.header.rows=oldRows;
  248.             } /*end if pairs<4*/
  249.         else { /*npair>=4*/
  250.             Spline( xd,yd, npair, c, x0,&y0, TRUE, am,bm,cm,dm,em,fm, &last );
  251.             for( i=1; i<=newRows; i++ ) {
  252.                 x0 = table.header.start + table.header.samp*((float)(i-1));
  253.                 Spline( xd,yd, npair, c, x0,&y0, FALSE, am,bm,cm,dm,em,fm, &last );
  254.                 table.header.rows = newRows;
  255.                 SetTable( i, col, y0, FALSE );
  256.                 table.header.rows = oldRows;
  257.                 } /*end for i*/
  258.             } /*end if pairs>=4*/
  259.     } /*end for col*/
  260.     
  261.     table.header.interpolated=TRUE;
  262.     table.header.rows = newRows;
  263.     Header2Vars();
  264.     CreateCol1();
  265.     DisposPtr(xd);
  266.     DisposPtr(yd);
  267.     DisposPtr(c);
  268.     DisposPtr(am);
  269.     DisposPtr(bm);
  270.     DisposPtr(dm);
  271.     DisposPtr(em);
  272.     DisposPtr(fm);
  273. }
  274.        
  275.  
  276.  
  277. TRIDI( A,B,C, D, T, N, E,F ) /* solves tri-diagonal matrix equation */
  278. float A[], B[], C[], D[], T[], E[], F[];
  279. int N;
  280.  
  281.     /*THE MATRIX EQUATION SOLVED IS :
  282.  
  283.     ( B1 A1 0  0  0  0  .  ) ( T1 )    ( D1 )
  284.     (                      ) (    )    (    )
  285.     ( C2 B2 A2 0  0  0  .  ) ( T2 )    ( D2 )
  286.     (                      ) (    )    (    )
  287.     ( 0  C3 B3 A3 0  0  .  ) ( T3 ) =  ( D3 )
  288.     (                      ) (    )    (    )
  289.     ( 0  0  C4 B4 A4 0  .  ) ( T4 )    ( D4 )
  290.     (                      ) (    )    (    )
  291.     ( .  .  .  .  .  .  .  ) ( .  )    ( .  )
  292.     (                      ) (    )    (    )
  293.     ( 0  0  0  0  0  CN BN ) ( TN )    ( DN )
  294.     
  295.     E, F are temp arrays of length N+2 */
  296.  
  297. {
  298.     int I, J;
  299.     float TEMP;
  300.         
  301.     E[1] = -A[1]/B[1];
  302.     F[1] =  D[1]/B[1];
  303.     for( I=2; I<=(N-1); I++ ) {
  304.         TEMP = B[I] + C[I] * E[I-1];
  305.         E[I] = - A[I] / TEMP;
  306.         F[I] = ( D[I] - C[I]*F[I-1] ) / TEMP;
  307.         } /*end for*/
  308.  
  309.     T[N] = (D[N]/C[N]-F[N-1]) / (E[N-1]+B[N]/C[N]);
  310.  
  311.     for( J=1; J<=N-1; J++ ) {
  312.         I=N-J;
  313.         T[I] = E[I] * T[I+1] + F[I];
  314.         } /* end for */
  315. }
  316.         
  317.